function [u,c,l]=UtilityFlow(NU_C,NU_L,PSI_L,zw,chat,TypeFlag)
%% Preliminaries
if isscalar(zw) && isvector(chat)
    zw      =   zw*ones(size(chat));
end
if nargin<=5
    TypeFlag    =   'SeparateUtility';
end
%% Solve the Optimal Utility Flow
switch TypeFlag
    case 'SeparateUtility'
        % Consumption
        if NU_L==0
            c       =   (zw/PSI_L).^(1/NU_C);
        else
            Temp    =   ( zw.^(1+1/NU_L) )*( PSI_L^(-1/NU_L) );
            TempInd =   Temp>1e-6;
            % Corresponding to the Unemployment
            c       =   chat;
            % Corresponding to the Normal Employed State
            if any(TempInd)
                temp    =   Temp(TempInd);
                temp_chat ...
                        =  chat(TempInd);
                TempFun =   @(c) c - temp.*( c.^(-NU_C/NU_L) ) - temp_chat;
                c0      =   temp.^( NU_L/(NU_C+NU_L) );
                lb      =   max(temp_chat,0);
                ub      =   lb + 2*c0;
                c(TempInd) ...
                        =   VecFun_Bisection(TempFun,lb,ub,1e-7);
                
            end
            
        end
        % Labor
        l       =   (c-chat)./zw;
        % Utility
        if NU_C==1
            u       =   log(c) - PSI_L*( l.^(1+NU_L)/(1+NU_L) );
        else
            u       =   c.^(1-NU_C)/(1-NU_C) ...
                        - PSI_L*( l.^(1+NU_L)/(1+NU_L) );
        end
    case 'GHH'
        l       =   ( zw/PSI_L ).^( 1/NU_L );
        c       =   chat+zw.*l;
        % Utility
        if NU_C==1
            u       =   log(c - PSI_L*( l.^(1+NU_L)/(1+NU_L) ));
        else
            u       =   ( c- PSI_L*( l.^(1+NU_L)/(1+NU_L) ) ) ...
                        .^(1-NU_C)/(1-NU_C);
        end  
end

